;;......................................................................................
;; NCL Graphics: 
;;......................................................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

begin

  plotfmt = "pdf"
  shortname = "AIE_T_2006_Tnudge" 

  ca = "AIE_1.9x2.5_EXP000"
  cb = "AIE_1.9x2.5_EXP001"
  cc = "AIE_1.9x2.5_EXP006"
  cd = "AIE_1.9x2.5_EXP008B"
  ce = "AIE_1.9x2.5_EXP008D"
  cf = "AIE_1.9x2.5_EXP005"

  na = ca
  nb = cb
  nc = cc
  nd = cd
  ne = ce
  nf = cf

  basepa = "/glade/scratch/kaizhang/run/run_"
  basepb = "/glade/scratch/kaizhang/run/run_"
  basepc = "/glade/scratch/kaizhang/run/run_"
  basepd = "/glade/scratch/kaizhang/run/run_"
  basepe = "/glade/scratch/kaizhang/run/run_"
  basepf = "/glade/scratch/kaizhang/run/run_"

  ma = ".cam2.h0.2006.ANN.nc" 
  mb = ".cam2.h0.2006.ANN.nc" 
  mc = ".cam2.h0.2006.ANN.nc" 
  md = ".cam2.h0.2006.ANN.nc" 
  me = ".cam2.h0.2006.ANN.nc" 
  mf = ".cam2.h0.2006.ANN.nc" 

  vname  = "dst_a1"
  vunits = "ug/m3"
  vscale = 1. ;;1e-3

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  pc = basepc+cc+"/"
  pd = basepd+cd+"/"
  pe = basepe+ce+"/"
  pf = basepf+cf+"/"
  fa = ca+ma
  fb = cb+mb
  fc = cc+mc
  fd = cd+md
  fe = ce+me
  ff = cf+mf
  casa  = (/"A: "+na, \
            "B: "+nb, \
            "C: "+nc, \
            "D: "+nd, \
            "E: "+ne, \
            "F: "+nf  /)
  caseM  = (/pa+fa, \
             pb+fb, \
             pc+fc, \
             pd+fd, \
             pe+fe, \
             pf+ff  /)
  case  = caseM

  Ncase = dimsizes(case)

  plotname = "figure_"+vname+"_zonal_"+shortname

  print(" ") 
  print(" ") 
  print("..................................................................") 
  print(" ") 
  print(" ")
  print(" varname   : "+plotname)
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1)) 
  print(" plot name : "+plotname)
  print(" ")
  print(" ")


;;......................................................................................
;; Open netCDF file and read in dimensions
;;......................................................................................

   filename=case(0)
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)

   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))

;;......................................................................................
;; open wks 
;;......................................................................................
  wks   = gsn_open_wks (plotfmt,plotname)             ; open workstation


;;......................................................................................
;; colormap 
;;......................................................................................

  colormap = "BlWhRe"
  colormap = "WhiteBlueGreenYellowRed"
  colormap = "BlueYellowRed"
  colormap = "BlueYellowRedWhite" 

  gsn_define_colormap(wks,colormap) 


;;......................................................................................
;; resources 
;;......................................................................................
  res1                      = True                ; plot mods desired
  res1@gsnFrame             = False               ; advance page
  res1@gsnDraw              = False                
  res1@gsnSpreadColors      = True                 ; use full colormap 

  if(colormap.eq."WhiteBlueGreenYellowRed") then 
    res1@gsnSpreadColorStart  = 30                    ; start at color 10
    res1@gsnSpreadColorEnd    = 235                    ; end at color 96
  end if 

  res1@cnFillOn             = True                 ; turn on color
  res1@cnLinesOn            = False                ; turn on/off contour lines
  ;res1@cnFillMode           = "AreaFill"           ; default "AreaFill" "RasterFill"
  ;res1@cnFillMode           = "RasterFill"           ; default "AreaFill" "RasterFill"
  res1@lbLabelBarOn          = True               ; turn on/off the labelbar
  res1@lbOrientation         = "Vertical"        ; labelbar Orientation "Vertical"
  res1@lbLabelsOn            = True               ; turn on/off the label of labelbar
  res1@lbBoxMinorExtentF     = 0.40                ; area allotted to each box of the LabelBar    
  res1@lbLabelFontHeightF    = 0.02                                               
  res1@lbLabelAutoStride     = False                ; optimal labels         
  res1@cnLineLabelsOn        = False                ; turn on line labels
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  res1@cnLevels             = (/0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20/)
  res1@cnLevelFlags         = False
  res1@tiXAxisFontHeightF    = 0.020               ; XAxis Font
  res1@tiYAxisFontHeightF    = 0.018               ; YAxis Font
 ;res1@tiYAxisString         = "Approximate Pressure(hPa)"
  res1@tiYAxisString         = "Pressure(hPa)"
  ;res1@tiXAxisString         = "Latitude"
  res1@vpXF                  = 0.12     ; x location left
  res1@vpYF                  = 0.88     ; y location  up
  res1@vpWidthF              = 0.65     ; width
  res1@vpHeightF             = 0.35     ; height
  res1@gsnStringFontHeightF  = 0.02
  res1@tiMainString          = " "
  res1@tiMainFontHeightF     = 0.03
  res1@tiMainFontThicknessF  = 0.35
  res1@trYReverse            = True
  res1@gsnYAxisIrregular2Linear = True

  res1@tmYRMode             = "Automatic"          ; turn off special labels on right axis

;;......................................................................................
;; read in variables & create plot
;;......................................................................................

  plot = new(6,graphic)
  
  do i=0,Ncase-1
     filename=case(i)
     fi=addfile(filename,"r")
    ;print(" read "+case(i))

     freqi = fi->FREQI 
     var   = fi->$vname$ 
     PS    = fi->PS

     T     = fi->T
     hyam  = fi->hyam
     hybm  = fi->hybm

;...................................................................................
; unit
;...................................................................................

     R = 8.3143

     pres = T

     nlev = dimsizes(T(0,:,0,0))

     do k = 0,nlev-1
        pres(:,k,:,:) = ( PS(:,:,:) / 100. ) * dble2flt(hybm(k)) + 1013.25 * dble2flt(hyam(k))
     end do

     var = var * pres * 100. * 29. / ( R * T * 1000. )    ; kg / m3
     var = var * 1.e+09

     fvar  = var
 
     P0mb =1000.
     hyam = fi->hyam
     hybm = fi->hybm

     ; type of interpolation: 1 = linear, 2 = log, 3 = loglog
     interp = 2
     ; is extrapolation desired if data is outside the range of PS
     extrap = False

     ; create an array of desired pressure levels:
    ;pnew = (/ 1000.,925.,850.,800.,750.,700.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100./)

     ;; use cam5 l30 levels, find and use the standard levels 

     pnew = (/ 100., 121., 143.,  \ 
               168., 200., 233., 274., \ 
               300., 379., 446., 500., \ 
               610., 700., 763., 821., \ 
               850., 887., 913., 936., \ 
               957., 976., 993., 1000. /) 

     pnew = pnew(::-1) 

     ; note, the 7th argument is not used, and so is set to 1.
     pvar = vinth2p(fvar,hyam,hybm,pnew,PS,interp,P0mb,1,extrap)

    ;printVarSummary(pvar) 

     var3d = dim_avg_Wrap(pvar(time|:,lev_p|:, lat|:, lon|:)) ;;dim_avg_Wrap(fi->$vname$(time|:,lev|:, lat|:, lon|:))

     var3d!0 = "time"
     var3d!1 = "lev_p"
     var3d&lev_p = pnew
     var3d&lev_p@units = "hPa"
     var3d!2 = "latitude"
     
     var2d = var3d(0,:,:) 
     var2d@_FillValue = -999

    ;printVarSummary(var2d) 
    ;print(var2d&lev_p@units) 

     var2d=var2d*vscale
     res1@gsnLeftString = casa(i)
     res1@gsnRightString= vunits                


    ;print(var2d) 

     var2d = where(var2d.le.0.,-999,var2d) 

     if(i.eq.0) then
        vara = var2d 
     end if
     if(i.eq.1) then
        varb =  var2d
     end if

     if(i.eq.2) then
        varc =  var2d
     end if

     if(i.eq.3) then
        vard =  var2d
     end if

     if(i.eq.4) then
        vare =  var2d
     end if

     if(i.eq.5) then
        varf =  var2d
     end if

     delete(var2d)
     delete(var3d) 
     delete(hyam) 
     delete(hybm) 
     delete(freqi) 
     delete(var) 

  end do 



  P0mb =1000.
  hyam = fi->hyam
  hybm = fi->hybm

  ;;vara = where(vara.le.0.,-999,vara) 
  ;;varb = where(varb.le.0.,-999,varb) 
  ;;varc = where(varc.le.0.,-999,varc) 
  ;;vard = where(vard.le.0.,-999,vard) 
  ;;vare = where(vare.le.0.,-999,vare) 
  ;;varf = where(varf.le.0.,-999,varf) 

  varb = varb - vara 
  varc = varc - vara 
  vard = vard - vara 
  vare = vare - vara 
  varf = varf - vara 


  res1@gsnSpreadColors      = True                 ; use full colormap 

  if(colormap.eq."WhiteBlueGreenYellowRed") then 
     res1@gsnSpreadColorStart  = 20                    ; start at color 10
     res1@gsnSpreadColorEnd    = 195                    ; end at color 96
  end if 

  res1@gsnLeftString = "CLIM"
  plot(0) = gsn_csm_pres_hgt(wks,vara(:,:), res1) 

  res1@gsnLeftString = "EXP02-EXP01"
  res1@gsnLeftString = "UVT6H-CLIM"
  delete(res1@cnLevels) 
  res1@cnLevels             = (/-5,-2,-1,-0.5,-0.2,-0.1,0,0.1,0.2,0.5,1,2,5/)

  plot(1) = gsn_csm_pres_hgt(wks,varb(:,:), res1) 
  res1@gsnLeftString = "EXP02-EXP01"
  res1@gsnLeftString = "UV6HT1D-CLIM"
  plot(2) = gsn_csm_pres_hgt(wks,varc(:,:), res1) 
  res1@gsnLeftString = "EXP03-EXP01"
  res1@gsnLeftString = "UV6HT4D-CLIM"
  plot(3) = gsn_csm_pres_hgt(wks,vard(:,:), res1) 
  res1@gsnLeftString = "EXP04-EXP01"
  res1@gsnLeftString = "UV6HT16D-CLIM"
  plot(4) = gsn_csm_pres_hgt(wks,vare(:,:), res1) 
  res1@gsnLeftString = "EXP05-EXP01"
  res1@gsnLeftString = "UV6H-CLIM"
  plot(5) = gsn_csm_pres_hgt(wks,varf(:,:), res1) 





  resP                 = True
  resP@gsnPanelYWhiteSpacePercent = 1  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 1  ; draw panel with white space added
  resP@txString            = vname
  resP@gsnPanelLabelBar    = False                ; add common colorbar
  resP@lbLabelFontHeightF  = 0.012               ; labels size
  resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 


  
  gsn_panel(wks,plot,(/3,2/),resP)
  
  print("ps file : "+vname)
  
end





